During some calculations on BH/aug-cc-pVTZ (in a 4e, 68o) active space, I noticed SCI energies going below the FCI energy by a small amount (around 1e-5 Ha).

Eventually, I tracked this down to a small discrepancy in how bits are masked in phase_single_det and phase_double_det.
The problematic calculation is the mask ~(1UL << (n + 1)) + 1. It generates the following pattern:

n, mask:
0  fffffffffffffffe
1  fffffffffffffffc
2  fffffffffffffff8
3  fffffffffffffff0
4  ffffffffffffffe0
...
60  e000000000000000
61  c000000000000000
62  8000000000000000
63  ffffffffffffffff

Here, one would have expected 0 for n=63. The result is that all bits in the determinant are counted up for the number of permutations, instead of no bits.

Below is a small example that shows the impact of this problem for the hydrogen molecule.
When below 64 orbitals, PyCI and PySCF energies match well, to 2e-11 Ha.
When above 64 orbitals, there was a small discrepancy of around 4e-8 Ha; this is reduced to around 4e-11 Ha after the proposed change.
If we permute the orbitals to reverse order (which moves an energetically important orbital into the problematic region), the discrepancy increases to 1e-3 Ha.

Script to produce the table above:

import pyscf
import pyci

xyz = "H 0. 0. 0.; H 0. 0. 1.10"
basis = "dyall-v4z"

mol = pyscf.M(atom=xyz, basis=basis)
nelec = (1, 1)
nelcas = sum(nelec)
ncas = mol.nao
print(f"{nelcas=} in {ncas=}")

mf = pyscf.scf.RHF(mol).run()
mf.conv_tol = 1e-12
mf.kernel()
cas = pyscf.mcscf.CASCI(mf, ncas, nelcas)
cas.kernel()

h1, ecore = cas.h1e_for_cas()
eri = pyscf.ao2mo.full(mol, mf.mo_coeff, aosym='1').reshape(ncas, ncas, ncas, ncas)
ham = pyci.hamiltonian(ecore, h1, eri.transpose(0,2,1,3))

wfn = pyci.fullci_wfn(ham.nbasis, *nelec)
wfn.add_all_dets()
op = pyci.sparse_op(ham, wfn)
e_vals, e_vecs = op.solve(tol=1e-15)
print(f"{e_vals[0]=:16.14f}")